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For  the  advection-diffusion  equation,  the  characteristic-Galerkin  formulations  are  obtained  by 
temporal  discretization  of  the  total  derivative.  These  formulations,  by  construction,  are  Eulerian- 
Lagrangian,  a-d  therefore  can  handle  time-dependent  domains  without  difficu  .y.  The  Galerkin /least- 
squares  space-time  formulation,  on  the  other  hand,  is  written  over  the  space-time  domain  of  a 
problem,  and  therefore  can  handle  time-dependent  domains  with  no  implementational  difficulty.  The 
purpose  of  this  paper  is  to  compare  these  two  formulations  based  on  error  estimates  and  numerical 
performance,  in  the  context  of  the  advection-diffusion  equation. 


0.  Introduction 

The  advection-diffusion  equation  arises  in  several  systems  in  fluid  mechanics.  It  is  also  a 
model  problem  for  the  Navier-Stokes  equations.  Time-dependent  domains  are  also  frequent, 
such  as  in  problems  involving  free  surfaces,  two  liquid  interfaces  and  drifting  objects. 

Although  the  advection-diffusion  equation  is  linear  and  much  simpler  than  the  Navier- 
Stokes  equations,  it  is  a  challenge  to  numerical  analysts  when  the  diffusion  coefficient  is  small. 
In  the  finite  element  framework,  there  are  several  methods  to  circumvent  the  difficulty.  To 
give  a  few  examples,  we  can  mention  the  streamline-upwind /Petrov-Galerkin,  Galerkin/ 
least-squares,  discontinuous  Galerkin,  and  the  Lagrangian-Eulerian  schemes  based  on  the 
characteristic  method  (see  [1]  for  a  review).  Time-dependent  domains  add  difficulty  to  the 
problem  because  one  has  to  deal  with  moving  meshes  and,  depending  on  the  formulation, 
projections  and  mesh  intersections. 

The  characteristic  methods  are  derived  from  the  analytical  solution  of  the  advection- 
dissipation  equation. 

Correspondence  to:  Dr.  O.  Pironneau.  Universite  Paris  6,  Analyse  Numerique,  T  55-65/5.  4  place  Jussieu  Paris 
75005  France 
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({> ,  +  uV(l)  +  acj)  =  f  on  Q  =  O  ^  ]0,  T[ , 

and  many  schemes,  such  as  those  proposed  by  Chorin  [2],  Bristeau  and  Dervieux  [3],  and  the 
particle  methods  are  based  on  this. 

In  1979,  adaptation  of  these  methods  to  the  finite  element  framework  was  found  con¬ 
currently  by  Benque  et  al.  [4],  Douglas  and  Russell  [5]  and  Pironneau  [6].  In  simple  words, 
the  idea  is  to  discretize  the  total  derivate  D<^/Dr  in  time  instead  of  <^,.  Thus 

0 ,  +  uV(t>  -  u  A<^>  =  /  on  ^  X  ]0,  T] 
could  be  discretized  in  time  by 

-  [0"'^'  -  -  u"(jc)k)]  -  u  =  f”  on  17  . 

This  is  the  basic  idea.  It  can  be  interpreted  as  a  splitting  method  (one  step  of  transport  +  one 
step  of  diffusion  [7])  because  it  looks  as  if  one  convects  (i>  first  and  then  applies  the  diffusion; 
but  the  error  analysis  shows  that  there  is  more  to  the  method  than  just  this.  In  particular,  the 
method  is  unconditionally  stable  and  in  some  cases  so  accurate  that,  for  example,  the  rotating 
hill  problem  can  be  solved  in  only  one  time  step,  if  desired.  As  pointed  out  by  Bermejo  [8], 
there  are  also  a  number  of  features  that  are  common  to  this  method  and  the  particle  method; 
it  is  in  fact  a  particle  in  cell  method  in  which  the  cut-off  function  is  the  hat  function  of  the 
finite  element  method,  and  for  which  a  projection  step  is  performed  every  time  step. 

The  Galerkin/ least-squares  space-time  (GLS/ST)  formulation  for  fixed  spatial  domains 
was  applied  by  Hughes  and  co-workers  and  Johnson  and  co-workers  (see  e.g.  [9-11]),  to  a 
large  class  of  fluid  dynamics  problems,  including  compressible  flows.  Tezduyar  et  al.  [12, 13] 
applied  the  GLS/ST  formulation  to  problems  with  time-dependent  domains,  and  provided 
several  numerical  examples  from  incompressible  fluid  dynamics,  including  those  involving  free 
surfaces,  liquid  drops,  two-liquid  flows,  and  flows  with  drifting  cylinders. 

In  the  GLS/ST  formulation,  at  each  time  step,  one  must  solve  the  problem  over  a 
space-time  slab.  For  time  dependent  domains  this  space-time  slab  becomes  non-cylindrical. 

In  this  paper,  we  first  describe  the  characteristic-Galerkin  and  GLS/ST  methods  and  give 
proofs  of  their  stability,  and  partial  proofs  of  the  error  estimates  when  the  spatial  domain  is  a 
function  of  time.  Then  we  evaluate  these  methods  on  a  test  problem  involving  advection- 
diffusion  of  an  exponential  hill  in  a  rotational  flow  field. 


1.  The  characteristic-Galerkin  method 

Consider  the  convection-dissipation  problem.  Find  <(>  such  that 


4) ,  +  uV(f>  +  a<f)=  f  in  Q  =  {x,  t:  X  E  J2(r),  t  E  ]0,  T[}  , 
<^(jc,  0)  =  <^)  V)  V;cE/2(0), 


(1) 

(2) 
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()  =  g{x, !)  VatET  (1)=^  {x:[u{x,t)-  u{x,t)\.nix,t)<0,x^d{}(t)} 

VrEjO,  r[,  (3) 

where  u,  /,  0“  g,  {},  F  are  given;  n  is  the  outer  normal  to  d/2(t)  and  i;(jc,  t)  is  the  velocity  of 
the  boundary  point  x  E  at  time  t.  As  usual,  f  is  time  and  x  E  R‘‘. 

PROPOSITION  1.  Assume  0(t)  is  bounded  with  Lipschitz  boundary  Fit).  Let  n  be  the 
normal.  If  th""  &  L^D),  g  S  C\dQ)  and  a,  u,  u  ,  E  L^(Q),  f  E  L'(Q),  then  system  (I),  (3) 
has  a  unique  solution  in  C‘^(/2) 

PROOF.  The  proof  is  constructive  and  will  be  used  also  to  justify  the  algorithm.  Let 
be  the  solution  of 

±  Y(.r\  =  f  ^  ’ 

dr  '  [u(A!’(t),  t)  otherwise ,  ''  ' 

with  the  boundary  condition 


^(0  =  ^  •  (5) 

If  u  is  the  velocity  of  a  fluid,  then  X  is  the  trajectory  of  the  particle  that  passes  at  position  x  at 
time  t.  Problem  (4),  (5)  has  a  unique  solution  because  u  is  Lipschitz.  In  (4),  (5),  x,  t  are 
parameters  so  X  is  written  X{x,  t;  r);  it  is  also  the  characteristic  of  the  hyperbolic  equation 
0).  Now  we  note  that 


^  <l>(X{x,  t;  t),  tI^,  =  <t>  ,{x,  t)  +  uV(i>{x,  t) 


(6) 


Thus  (1)  is  rewritten  as 


(f>r  +  a<f)=f 


(7) 


and  it  can  be  integrated  analytically; 


(f>(x,t)  =  exp  a(X(T),T)dT  A  +  /(X((r),  a)  exp  a(X(T),  t)  dr  da  . 


(8)  - 


To  find  A  we  use  (2)  and  (3).  If  X^x,  t;  0)  E  F,  then 
A  =  g(X(x,  /;  0),  t(a:))  , 

where  t(a:)  is  the  largest  time  (<r)  such  that  Ar(A:,  t;  t)  E  F{t). 
If  X{x,  t;  0)  E  /2(0),  then  JDXin  ^ 


(9) 


ft-l  ic^  j 
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\  =  ti>\X{x,r,0)).  (9') 

Let  be  a  bounded  open  set  big  enough  to  contain  fl(t)  at  all  times.  Note  that 

V •  =  <pV •  u  +  uV<i> .  (10) 

Thus,  we  also  have  an  existence  theorem  for  the  convection-diffusion  equation  in  divergence 
form.  (Recall  that  //(div,  f2)  =  {uGL’(/2)":  VwGL^( /}))).  We  also  introduce  the  notation 

X  =  {x,  r;  X  G  r{t),  t  G  ]0,  T[}  . 

Let  us  assume  that  u  can  be  extended  into  such  that  u  G  L^(0,  T:  //(div. 

COROLLARY.  With  the  hypotheses  of  Proposition  1  and  if  u  G  L‘(0.  7; //(div,  Z2“))  n 
L"(Q),  problem  (/')-(3), 

<f>,-¥V-{u<j))=f  inQ  =  nx]0,T[,  (1') 

has  a  unique  solution  in  C‘^(0,  T;  L '(/})). 

PROOF.  Use  (10)  and  Proposition  1  with  a-V  u. 

We  also  recall  a  similar  result  for  the  convection-diffusion  equation. 

PROPOSITION  2.  Let  fl  be  bounded  with  boundary  F  Lipschitz.  The  problem 

, -1- V- (u<f>)  -  r' A<f>  =/  inQ,  (11) 

0(x,O)  =  </»V)  in  12(0),  (12) 

4>  =  g  on  2  ,  (13) 

has  a  unique  solution  in  L^  in  t,  in  x  if 

m.,u,^GL’=(Q),  feL^iQ),  cl>^^L\n),  gisL^int,  in  x , 

v>0.  (14) 

PROOF.  The  proof  can  be  found  in  Ladyzhenskaya  et  al.  [14]  and  Lions  [15]. 

Extensions  of  this  result  to  the  time  dependent  domain  can  be  made  by  using  conformal 
mappings  for  example  which  map  fl{t)  into  /2(0).  The  problem  then  is  on  a  fixed  domain  but 
the  coefficient  v  is  now  time  and  space  dependent. 

l.l.  Discretization  in  time 

Now  for  clarity  we  assume  that  V  -  n  =0.  By  (6), 
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<i>,  +  uVti>=~  MX{x,  /;  t),  t)U,  .  (15) 

Thus  by  using  the  fact  that  X(x.  (n  +  l)k-,{n  +  l)/c)  =  .v.  we  can  say  that 

^  [r^\x)  -  nx^ix))] ,  (i6) 

where  X''(x)  is  an  approximation  of  A!'(x,  (n  +  \)k-,  nk). 

Let  us  denote  by  X"  an  O(k^)  approximation  of  and  by  A'"  an  approximation  of  order 
0{k^)  (the  differences  between  the  indices  of  X  and  the  precision  order  are  due  to  the  fact 
that  X"  is  an  approximation  of  A'  on  a  time  interval  of  order  k;  a  scheme  of  Oik")  gives  a 
precision  of  0(A:“^')). 

For  instance. 


Xlix)  =  x-u''ix)k,  (17) 

X^ix)  =  X  -  u''*^''i^x  -  u’'(x)^^k  ,  (18) 

modified  near  the  boundary  without  losing  accuracy  (see  below)  but  so  that  A';(r2"*')Cf2''. 
Then  we  derive  a  scheme  for  (11): 

^  (0"^' -  <^"0,^")- .  (19) 

Note  that  0"^’  and  <f>''°X"  are  both  defined  on  n”*'. 

REMARK  1.  A  second  order  scheme  could  be 

i  -  d^^oxl)  -  ^  A(<f.''"‘  +  V)  =  ,  (20) 

but  there  is  a  difficulty  because  and  eft"  are  not  defined  on  the  same  domain. 


We  will  not  prove  convergence  in  the  general  case.  For  simplicity,  we  will  assume  that  the 
distance  between  H”  and  A'i(f7"^‘),  8(0",  Xlifi”*'))  is  0(A:^).  This  is  not  possible  unless  the 
normal  velocity  of  the  fluid  is  equal  to  the  normal  velocity  of  the  boundary, 

u  -  n  =  V  n  on  2  • 


LEMMA  1.  Ifuis  regular  and  if  Xl{fl''^')C.  the  schemes  (19)  and  (20)  are  L^-siable  and 
if  Sin",  XKn''^'))  =  O(k^),  U  is  convergent  0(k). 


PROOF.  We  multiply  (19)  by  and  integrate  over  Unless  explicitly  mentioned  all 

norms  and  integrals  are  on 


n  +  1 1 


(21) 
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Now  by  a  change  of  variable  y  =  AT'JCa:),  we  find  that 

=  det|vx:r'  dy  «  (22) 

Therefore  0"  satisfies 

l|0"IL«c[|/|„,;,  +  |<j"Ul.  (23) 

To  obtain  the  error  estimate,  we  proceed  exactly  in  the  same  way  but  on  the  difference 
e"  =  0"  -  instead.  Indeed 

i  («"*■  - c^xi)-vS€-'  =  I  ^  (("")  =  o(*:) , 

when  u  •  n  =  u  •  n  on  the  boundary.  Thus,  (23)  gives 

lk;i|sc[IU“IL.n»-lo(*)ll. 

1.2.  Approximation  in  space 

Let  us  now  discretize  (19)  or  (20)  with  the  finite  element  method.  Then  we  obtain  a  family 
of  methods  for  which  no  additional  upwinding  is  necessary;  the  schemes  are  unconditionally 
stable! 

For  example  with  (19),  assuming  that  u  has  zero  divergence,  a  possible  scheme  would  be 

•(-r 'H-r  ■  + )cr  vrrv^r  dx = ,t:(x:(x))H^r(x)  dx 

+t dx  (2“) 

where  //qa  *  is  a  space  of  continuous  polynomial  functions  on  a  triangulation  of  with 
zero  traces  on  the  boundaries  and  is  a  polynomial  approximation  inside  fl  of  g. 

Notice  that  (24)  is  a  linear  system  of  the  type 

=  Agi,  +  b  ,  where  /t,,i  "'V'  dr: . 

b=\^„{kf  +  r,mx))\«Xx)dx  , 
and  where  {w‘}  is  a  basis  of  //q/T'  ^nd  =  E 

REMARK  2.  A  change  of  variable  can  be  made  in  the  last  integral  of  (24),  x-*  y  =  K"(x),  so 
(24)  is  equivalent  to 
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dx  +  dJt = k  /-'.v,,  (lx 

+  (  My)>^A(-^':r'(y))dciVX"dy  Vh-, e *r' e H*" 


(24') 

PROPOSITION  3.  If  V-  «  =0,  u  •  All,.  =  V  ■  «lr,  scheme  (24)  is  L\n)  stable 

even  if  v  =  0. 

PROOF.  Replacing  w,,  by  in  (24), 

liar'll: i^r'pdx +  dx 

=  (ij:+ .«.:(.v:(x)),(.r'(x)dx 


(khlo(l  +  clulf.,,^)  +  /c|/""‘!„)l(/)"^" 


(25) 


The  last  inequality  is  a  consequence  of  the  fact  that  the  mapping  x -»•  A'  is  volume  preserving. 
Finally, 


(26) 


l<A:io,n  ^  (1  +  c|«|, +  S 


PROPOSITION  4.  If  Hll'  is  the  space  of  continuous  piecewise  affine  functions  on  the 
triangulation  ofQ"*',  then  the  L'iO)  error  between  </>"  solution  of  {24)  and  (f"  solution  of  {II ) 
is  0(h^lk  +  h).  Hence  the  scheme  is  0{h'lk  +  k  +  h). 


(27) 


PROOF.  Let  us  subtract  (11)  from  (24)  to  obtain  an  error  projected  on  L‘: 

rt+l  1^4-1  f-T  j  n  +  l 

=<Ph  ~  ^h<P 

where  is  an  interpolation  in  of  We  obtain 

dx + k.  £...  v.r'vw,  dx  -  .i-xy,  dx 

=  £..,(d>"*‘  -/7,(<>"*')»'»dx+  V((i)"*'  -  n,(l)"*')VB.,dx 

~  f. -n>r)-x’;w,dx.  (28) 


From  (28),  with  w,,  =  we  find 
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<(i  +  c|m|i  ,A')||€^I1^  +  ( I  +  ' 


n  +  I 

h  V  ^ 


SO 


lur'il.  «(1  +  c|«l?..*:)IU:lL  +  Qh^  +  M) .  (29) 

REMARK  4.  By  comparing  with  (24),  we  see  that  is  solution  of  a  problem  similar  to  that 
of  <t>^  in  which  /  would  be  replaced  by 


-n,r 


')-^{r-n,r)ox:. 


where  A,,  is  an  approximation  of  A.  We  have  bounded  the  first  and  the  last  terms  in¬ 
dependently.  By  being  more  cautious  one  can  show  that  the  error  is  in  fact  0{h  +  k  +  min(fi7 
k,h))  [5]  but  the  dependence  on  v  is  stronger. 

The  case  v  =  0:  Then  (24)  becomes 


=  r.oX'lw.dx  +  k  dx 

vw,  <^r‘e 


(30) 


that  is, 


K'-nMl^xi  +  kr^),  (31) 

where  77^  is  the  projection  operator  from  into  From  this  equation  we  see  that  the 

fewer  time  steps  there  are  the  better  because  a  projection  inevitably  produces  numerical 
diffusion;  however  when  k  is  large  the  computation  of  X'’  is  costly.  Experience  shows  that 
k=^l.5hlu  is  a  good  choice.  This  is  a  striking  behavior  of  these  Lagrangian / Eulerian 
methods:  they  should  be  used  wun  a  CFL>  I. 

Conservativity.  Assume  that 

u  -  n  =  V  ■  n  on  2  ,  (32) 

and  that  /  =  0  and  v  =  0;  the  equation  is 

4>  ,  +  uV4>  =  0,  (t>ix,  0)  =  <t>\x) .  (33) 

By  hypothesis  V  -  «  =  0,  so  by  integrating  (33), 

I  <f>(f,  x)dx=  [  (t>\x)  dx  Vr .  (34) 

Jou)  Jsi(O)  ^  ' 

On  the  other  hand,  (24')  with  =  I  also  gives  (34)  because  det|VA"’|  =  I  and  A"'(/2''^')  = 
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12".  So  the  scheme  is  conservative  up  to  quadrature  errors  on  the  approximation  ot  the  last 
integral  in  (24)  or  (24’)  and  up  to  integration  errors  due  to  X"  being  replaced  by  A'"  or  A"!. 
Indeed  W;,  =  I  in  (24)  gives 


dj:  = 


dt^^X"  dx  = 


«^:(y)detlVr;r'dy  . 


(35) 


Thus,  if  det|VAr"|  =  1  (which  implies  A'"(/}''"‘)  =  /2"),  we  obtain 


.  ‘  (^)  dJC  =  d>:(y)dy  =  ^'V)  d^:  . 


(36) 


Thus,  the  scheme  is  conservative  up  to  numerical  integration  errors  on  X".  In  many 
applications,  this  may  not  be  sufficient  but  (17)  or  (18)  can  be  modified  so  that  det|TA'"|  =  1. 
Indeed  since  V-  u  =  0.  there  exists  a  stream  function  (vector  in  3D)  <//  such  that  u  =  V  x  (//.  Let 
be  a  continuous  piecewise  approximation  of  l^.  Then  the  ODE  dA'/dT  =  V  x  can  be 
integrated  exactly  ({A’(t)}  is  a  broken  line)  and  jc— ►  A'"(x)  is  volume  preserving. 


1.3.  Implementation  problems 

Two  points  remain  to  be  clarified:  (i)  the  computation  of  Ar''(.v)  and  (ii)  the  computation  of 
the  integral, 


=  (37) 

1.3.1.  Computation  of  (37) 

Primal  Gauss  formula.  We  use  a  Gauss  quadrature  formula, 

(38) 


for  instance  with  elements  we  can  take: 

(a)  =  all  the  mid-points  of  the  edges  and  w*  =  cr^/3  in  2D,  o-^/4  in  3D,  where  is  the 
area  (volume)  of  the  elements  which  contain 

(b)  the  3  (4  in  3D)  inner  nodes  Gauss  formula  [16, 17]; 

(c)  the  7  integration  points  formula  (vertices,  mid-edges  and  central  node). 

Any  other  integration  formula  can  be  used  as  long  as  the  weights  are  positive.  Numerical 
tests  can  be  found  below  with  triangles  and  in  [18]  with  quadrangles  and  the  9  nodes  formula. 

Dual  Gauss  formula.  Denote  by  X"^*  {y)  the  inverse  mapping  of  =  A'^(jr).  Obviously 
A'^^(y)  is  an  approximation  of  the  solution  at  (n  +  \)k  of  dAT/dr  =  u(A'(t),  t),  X{nk)  =  y. 
Another  class  of  methods  can  be  obtained  if  a  change  of  variable  is  made  in  the  integral  as  in 
(24'). 


djc 


.  j 


(39) 
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When  V  -  u  =  0.  the  term  det|VA'^,|  is  1.  so  a  Gauss  formula  gives 

(41)) 

Now  with  such  a  quad;:  .are  formula,  the  scheme  is  conservative  regardless  of  quadrature 
errors.  Such  schemes  were  introduced  and  thoroughly  tested  by  Benque  et  al.  [4).  However, 
the  effects  of  quadrature  errors  on  stability  and  error  estimates  is  an  open  problem. 

1.3.2.  C  ('  npuiation  of  X^ix) 

To  compute  one  must  answer  the  following;  given  find  /  such  that 

r, ,  a  triangle  of  the  triangulation  of  D".  This  is  not  a  simple  problem.  We  shall 
proceed  in  two  steps. 

Step  1.  Find  m  such  that  E  T^,  a  triangle  of  the  triangulation  of  fi".  In  our  problem  we 
know  only  in  which  triangle  of  /}'’*'  lies;  we  do  not  know  that  information  in  the 
triangulation  of  f7".  So  assume  we  know  {m'.  with  f  *  near  and  E  T,„  ,  a 
triangle  of  the  triangulation  of  fl".  Then  intersect  the  segment  {f'\  f*}  with  the 
triangulation  of  ii"  and  proceed  from  neighbor  triangle  to  neighbor  triangle,  starting 
from  4^'*  until  is  reached.  This  is  possible  if  we  have  an  array  giving  the  number  of 
the  three  neighbor  triangles  of  each  triangle  of  the  triangulation  of  iT. 

Step  2.  Compute  all  intersections  of  with  the  edges  of  the  triangulation  and 

proceed  from  neighbor  to  neighbor  until  is  reached  to  find  /, 

But  then  at  no  extra  cost  we  can  improve  on  the  scheme  to  compute  A"’(  ):  the  numerical 

scheme  of  X''(^'‘)  is  applied  within  an  element  and  when  another  element  is  reached,  the 
value  of  on  the  new  element  is  used.  With  the  first  order  scheme  then,  one  can  use  the 
barycentric  coordinates  {A,}  of  and  proceed  as  shown  in  Fig.  1.  First  on  each  element 
compute  fx^  such  that 

u(Q)=  2 

/  =  1 . 

where  {q‘}  are  the  vertices  of  the  element,  Q  its  barycenter  and  d  the  dimension  of  the  space. 

Then  find  p  such  that 


Fig.  1.  Computation  of  XI  by  following  the  streamline  in  the  triangulation. 


().  Fironneau  et  ciL.  C  haractenstic-iialerkin  and  italerkini  k'usi-squures  \i)ai.i’-iirnt'  !ornnauiii)n.s 


^  A'Sst).  (42) 

To  tind  p,  we  assume  that  m  is  the  index  ot  the  A,„,  which  is  zero; 

P=-—  ,  (43) 

F-m 

then  if  A,  ^0,  V/.  the  m  is  the  right  one. 

Thus  most  of  the  work  goes  into  computing  p.,  for  all  the  elements.  Because  of  round  off 
errors,  it  may  be  difficult  at  times  to  decide  which  is  the  next  triangle  that  the  characteristic 
will  cross,  for  instance  if  is  on  a  vertex;  so  careful  programming  is  needed. 

If  =  V  X  \ii^  and  il/t,  is  piecewise  affine  and  continuous,  as  mentioned  above,  a  slightly 
different  algorithm  can  be  used  in  place  of  (17).  Instead  of  defining  the  charactenstics  as  the 
curves  tangent  to  u,, ,  we  define  them  as  curves  of  equal  values  of  ;  the  numerical  advantage 
is  that  if  a  characteristic  enters  a  triangle  if  must  leave  it  even  in  presence  of  round  off  errors, 
i.e,  given  an  entry  point  and  one  needs  just  to  test  which  side  q'q'  is.  such  that 

ilj,,{q  )E.[4/,J^q‘u),  ilt^,{q’)\  or  <A,.( p' )]  and  then  obtain  the  exit  point  by  linear 

interpolation  based  on  0,,. 


2.  Galerkin  least-squares  space-time  formulation 

Consider  again  the  convection-diffusion  equation 


(t) ,  +  V-  -  V  =  / 

on  Q  , 

(44) 

<f>(x,0)  =  <t)^\x) 

in  42(0) , 

(45) 

=  8 

on  2  • 

(46) 

Assume  for  simplicity  that  V  •  «  =  0. 

Choose  a  time  step  k  and  construct  a  quadrangulation  Tl  of  /if/"),  t"  =  nk,  for  each  n  such 
that  the  quadrangulation  of  is  obtained  by  moving  the  vertices  of  that  of  f2(r") 

with  the  condition  that  no  quadrangle  is  flipped  over  in  the  process.  This  way,  by  joining  the 
vertices  of  with  their  images  in  one  obtains  a  proper  quadrangulation  of  the  space 

time  slab 

Q""'  =  {{X,  ty.xen{t),te]t\  r""‘[}  .  (47) 

REMARK.  The  method  does  not  require  that  the  mesh  at  the  next  time  step  should  be 
obtained  by  moving  the  vertices  of  the  previous  mesh,  but  our  proof  of  convergence  relies  on 
this  special  case. 

Define  the  space  of  piecewise  linear  functions  which  are  piecewise  linear  in  each  component 
of  the  vector  (x,  t)  of  which  are  continuous  in  x  but  discontinuous  in  r: 
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C  (Q""'  ):  piecewise  affine  in  x^,  affine  in  /  on  Ql)  . 


W"'"'  = 

Oh 


[w.EW;;"';  w,=0on  E}  • 


(48) 


The  Galerkin  least-squares  space-time  (GLS/ST)  method  defines  the  approximate  solution 
as  the  solution  in  equal  to  on  E,  and  such  that  for  all  in  we  have 


t.=  I  '  y  ’ 

(w.  , +  u-Vw^-  i^AwJ  +  (<^r' 


where  n^,  is  the  number  of  elements.  For  the  purpose  of  proving  well-posedness  and 
convergence,  we  use  the  following  modified  form  of  (49): 

Ig..,  !'<>».:'  +  +  "^‘^*>1 +/„,  ('*r'  - 

+  ‘'/g..,^<''r'V^,=0.  (50) 

The  spatial  domain  i}”  is  to  be  understood  as  Z}”  x  (r"}.  The  scalar  t  is  a  positive  parameter 
of  order  0{h  +  k). 


PROPOSITION  5.  Problem  {50)  is  well-posed',  the  solution  exists  and  is  unique. 


PROOF.  Equation  (50)  is  a  linear  system  with  respect  to  the  values  of  at  the  vertices  of 
the  quadrangulation  of  Q"*',  i.e.  the  vertices  of  both  Tl  and  T^^'.  There  are  obviously  as 
many  equations  as  unknowns  so  we  only  need  to  check  that  the  kernel  of  the  system  is  zero, 
i.e.,  that  f=g  =  4>^  =  0  implies  =  0  for  all  n. 

With  these  zero  data  and  w^  =  ^  (50)  is 


L-  [5  +  \  +  /„.  wr')'  -  KK'] 

+  w  T = 0 + '  + u  .v,fr '  )^ 

Let  us  integrate  the  first  two  terms  by  parts; 

')^i = (K'f  - 1.  (■*:*')= 

■"I-  Lo 


(51) 


(52) 
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Now  the  last  integral  is  positive  because  is  zero  when  u  ■  n  <  v  ■  n.  Therefore  (50)  yields 


Define 


(53) 

(54) 


Then,  using  Schwartz’  inequality,  (53)  is  also 

But  this  can  be  rewritten  as 


(55) 


therefore 


Now  =  0  by  hypothesis  and  a„  5^0,  so  (56)  implies  a„  =  0. 


(56) 

(57) 


Let  us  prove  convergence.  For  simplicity  we  will  do  it  in  the  hyperbolic  case  only.  When 
v^Q,  the  method  converges  0(h  +  k)  in  H\  Extension  to  the  general  case  can  be  established 
following  [11]. 

PROPOSITION  6.  When  r  =  0(/i  +  ^),  v  =  0,  g  =  the  method  converges  0{{k  +  h)^'^)  in 

l\ 


PROOF.  The  proof  is  an  adaptation  of  that  given  by  Johnson  [11].  The  idea  is  not  to  treat 
space  and  time  separately.  Let  us  introduce  the  vector 


u  =  {i,uyeR'‘^' 

(58) 

and  rewrite  (44),  for  =  0,  as 

, <!>=/■ 

(59) 

Thus  (50)  is  equivalent  to 

-/![»'» +  [rr-  = o . 

(60) 

The  only  difference  with  the  time  independent  formulation  of  the  least-square  Galerkin 
method  is  the  fact  that  the  elements  are  discontinuous  in  time  and  the  presence  of  the  last 
integral  in  (60)  which  looks  like  a  penalty  term  to  insure  that  is  close  to  at  t  =  t”. 
Therefore,  let  /7^  be  an  interpolation  operator  from  H\Q)  into  11  Wl  and  define 
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=  n^<t>  ,  e  *  -  <t>^  . 


Subtract  (59)  in  variational  form  from  (60): 


I 


[w^  +  rUV,_,w^ 


(61) 


(62) 


Notice  that 


<*«/,-  ('/'a  -  • 

Choose  Wf,  =  el*'  and  rewrite  (62)  as 


Jj..,  (uv,., 

-  f^..,  -  r-')ur'  *ruv„,t:-'] 

+l.m"-r-’)-w-r)ur^ 

Next  notice  that  if  N  denotes  the  space-time  normal  to  BQ,  we  have 


(63) 


(64) 

(65) 

(66) 


The  last  integral  is  positive  because  e^  is  zero  when  w  n<v  n.  Using  this  in  (63)  and 
integrating  by  parts  the  term  UV^  X'I'l*'  -  ,  yields 


(67) 


Finally  we  rewrite  the  first  three  terms  as 
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jn  JSI  Jil 


and  we  obtain 


dSiU) 


jkr'ILr— |KlLr  + 1 1  -  ^-'r  +  H- 

+  |r(f/V,,er')'U- 

S  («”  -  +  £.  tK  -  -  n*') 

+r  f  (K"  '^‘V 

Ji''  JdOd) 

+  CT|f/V^  X'Ar' ~ ''  r"  ■ 


-  v)-n 


(69) 


Here  we  have  used  the  fact  that  =  tl/"  and  on  /2"  (see  (61)).  It  would  not  be 

true  if  the  triangulation  at  step  n  +  1  was  not  obtained  by  moving  the  nodes  of  that  at  step  n. 

Now  denote  =  {(jc,  t):  x  G  dn(t),  t  ^  r)  and  =  {(x,  r)  G  0:  r «  T}.  Let  us  sum  with 
respect  to  n: 


n=0 

= w:\in« + 21*.  -  + 2  £„ « -  *•”).: 

-  2  +  2  2  l*/'h  “  lo.n" 

+  2  [  (lA;,  -  </))€;,(«  -  v)-  n  +2T\UV^  ,(h  ~ 


(70) 


Now.  since  t/r,,  is  an  interpolate  of  d».  all  terms  on  the  right  are  bounded  by  some  power  of 

.  3/2 


8  ~  h  +  k,  the  biggest  of  which  is  5  because 

\^h  -  <^'lo.s^c||(/)||2_e5^'^ 

on  any  smooth  surface  5  of  So  the  right-hand  side  of  (69)  is  bounded  by 


(71) 

(72) 
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+  5(2;  [e:  -  .r'lr, ' 

*-  '  /I 

+  -  u)- ne™!^  j-  .  (73) 

Assume  t  =  0(5 ).  Since  =  0  when  (i/  -  y)  •  n  <  0.  we  have 

2  k"  -  +  f  e\u-v)-n  +  T\UV^,€\l^^ 

<c|l0||,.g5^'^  (74) 

REMARK.  If  II  0(1 2  is  replaced  by  ((0II3,  then  5^^^  can  be  replaced  by  5'  in  (72).  With  this 
additional  regularity  and  our  hypothesis  the  method  is  0(5 “). 


3.  Numerical  simulations 


3.1.  Summary  of  the  algorithms 

The  characteristic-Galerkin  method  used  is  based  on  a  finite  element  method  of  order  1 
(continuous  piecewise  affine)  on  triangles  with  a  primal  Gauss  quadrature  with  seven 
integration  points,  the  vertices,  the  middle  of  the  edges  of  the  triangles  and  the  center  of  the 
triangles.  Finally  the  characteristics  are  computed  exactly  via  a  piecewise  linear  stream 
function. 

Since  the  domain  is  time  dependent,  the  triangulations  at  each  time  step  are  different.  Here 
we  assume  that  each  triangulation  is  obtained  from  the  previous  time  step  by  moving  the 
vertices  with  the  velocity  field  v.  We  have  two  velocity  fields,  the  mesh  velocity  v  and  the  fluid 
velocity  u.  At  each  time  step  the  following  must  be  done; 

ALGORITHM  1 :  Characteristic-Galerkin  with  time  dependent  domains 

(1)  Compute  the  coordinates  of  the  vertices  of  the  new  mesh. 

(2)  Compute  the  old  mesh  triangle  numbers  which  contain  the  new  mesh  quadrature  points. 

(3)  Find  where  these  quadrature  points  were  in  the  fluid  at  the  previous  time  step  and  their 
triangle  number  in  the  old  mesh. 

(4)  Compute  the  integrals  just  on  the  right  of  the  equal  sign  in  (24)  with  w,,  =  (new  mesh); 
use  the  seven  Gauss  point  formula. 

(5)  Compute  the  matrix  of  the  linear  system  (>v',  w'  are  with  respect  to  the  new  mesh) 


i  i  \ 

w  w  + 


Vvv'Vw' . 


(75) 


(6)  Solve  the  linear  system  (by  the  Choleski  algo  ithm  for  instance). 

(7)  Update  all  arrays  and  go  back  to  1  until  /  =  T  is  reached. 
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It  should  be  noted  that  when  =  0  and  when  the  motion  of  the  mesh  is  incompressible 
iy  ■  V  —  0),  the  matrix  ^4  does  not  depend  on  n.  This  allows  us  to  program  the  method  to  have 
always  only  one  mesh  at  a  time  in  the  memory.  It  also  allows  a  unique  factorization  of  A  but 
we  have  preferred  to  reconstruct  A  at  every  time  step. 

REMARK.  Algorithm  1  computes  the  solution  of  (11)  in  the  frame  of  reference  in  which  the 
domain  is  fixed  so  that  the  velocity  of  the  mesh  should  be  changed  accordingly  if  the  results 
are  desired  in  the  original  frame  of  reference. 

The  least-squares  Galerkin  method  is  conceptually  simpler  and  much  easier  to  program,  but 
it  requires  the  solution  of  a  linear  system  which  is  double  the  size  of  the  other  method  and  the 
system  is  non-symmetric.  On  the  other  hand  it  is  exactly  conservative  (the  other  method  is 
conservative  only  if  dual  quadratures  are  used). 

ALGORITHM  2:  Least-squares  Galerkin 

(1)  Triangulate  Q". 

(2)  Construct  the  linear  system  of  matrix 

A,  i  =  +  iNw‘)(w'  +  +  uVw')  ■+•  vVw'Vw^]  -i-  ^  w'w'  .  (76) 

(3)  Construct  the  right-hand  side  and  solve  the  system. 

Contrary  to  Algorithm  1.  the  results  of  Algorithm  2  correspond  to  the  solution  of  the 
problem  in  the  fixed  frame  of  reference  where  the  mesh  moves. 

3.2.  The  rotating  hill 

The  flow  is  purely  rotational  and  the  convected  variable  has  an  exponential  hill  profile,  i.e.. 
almost  zero  everywhere  excent  in  a  small  region.  After  one  turn,  the  solution  should  be 
identical  to  the  initial  condition,  if  j/  =  0. 

The  velocity  of  convection  is 


u  =  (o  X  X  =  {10X2,  -  wx, )'  ,  a>  =  1  . 


The  initial  condition  is 

0“(jc)  =  exp(-lO|(x-.r„)h,  =  (0.5  0.0)'. 

The  domain  Q  is  the  unit  circle.  The  time  step  chosen  is  6/  =  -rr/M. 

Experiments  are  done  with  a  moving  domain  and  the  mesh  moves  with  the  velocity  of  the 
domain  v  which  is  different  from  the  *  '.’cif'  of  the  fluid  u. 

u  =  (u,,  Uj)' ,  u  =  P{x,t)xx, 

i8(jc,  t)  =  I3„  +  co&(2'nt)(3^  sin(2TTar) ,  r  =  |j:|  . 
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Honnini*!  Aais 


V«<ite«l  Axis 


Exact  aotutioA 


Compulabon  don*  on  •  fix*d  domain 


Honiontai  Asia 

Computation  with  a  moving  domain 

Fig.  2a.  Central  cross  sections  of  the  exponential  hill  with  «/  =  0  obtained  with  the  characteristic-Galerkin  method 
after  14  time  steps;  exact  solution  (top),  for  fixed  mesh  (middle),  and  for  rotating  mesh  (bottom). 
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Fig.  2b.  Contour  plots  of  the  exponential  hill  with  =0  obtained  with  the  characteristic-Gaierkin  method  after  14 
time  steps;  initial  condition  and  mesh  (top),  for  fixed  mesh  (middle),  and  for  rotating  mesh  (bottom). 


At  each  iteration  in  time,  the  vertices  q'  of  the  triangulation  are  moved  by  integrating  the 
equation 


=  y (x(f),  t) ,  xit°)  =  q  , 


with  ten  steps  of  the  forward  Euler  scheme. 
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Fig.  2*.  Contour  plots  of  the  exponential  hill  with  v  =  0  obtained  with  the  least-square  Galerkin  method:  initial 
condition  (top),  for  fixed  mesh  (middle),  and  for  rotating  mesh  (bottom). 


We  tested  several  cases  in  which  values  of  the  parameters  in  the  definition  of  )3(jc,  t)  are 
^0  =  0  or  -1,  =  0  or  0.25,  a  =  0  or  1  and  v  =  0  or  0.01.  In  all  cases,  the  solutions  obtained 

with  the  characteristic-Galerkin  and  the  GLS/ST  methods  are  presented.  We  use  a  time  step 
size  of  rr/ 14. 

Rotating  hill  with  fixed  and  rotating  meshes.  We  set  =  0  for  the  fixed  mesh  and  /3„  =  - 1  for 
the  rotating  mesh.  We  computed  in  each  case  for  14  time  steps.  Figures  2a,  2b  and  2*  show  the 
solutions  for  the  case  with  i'  =  0  at  time  steps  7  and  14.  while  Figs.  3  and  3*  show  the  solutions 
for  the  case  with  v  =  0.01  at  the  same  time  steps. 
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Fig.  3.  Central  cross  sections  and  contour  plots  of  the  exponential  hill  with  t'  =  0  01  obtained  with  the 
characteristic-Galerkin  method  after  14  time  steps:  for  fixed  mesh  (top),  and  for  rotating  mesh  (bottom). 


Fig.  3*.  Contour  plots  of  the  exponential  hill  with  v  =0.01  obtained  with  the  least-square  Galerkin  method;  for 
fixed  mesh  (top),  and  for  rotating  mesh  (bottom). 
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Rotating  hill  with  stretching  mesh.  In  this  case,  we  set  13,  =  0.25  to  describe  the  motion  ol  the 
mesh.  The  viscosity  is  set  to  U.  Figures  4  and  4*  show  the  meshes  and  the  solutions  at  time 
steps  7  and  14. 

Rotating  hill  with  rotating  and  stretching  rnesh.  In  this  case,  we  set  (3„  -  -  1  and  /j,  =  0  25  to 
describe  the  motion  of  the  mesh.  The  viscosity  is  set  to  0.  Figures  5  and  5*  show  the  meshes 
and  the  solutions  at  time  steps  7  and  14. 

The  computations  for  the  characteristic-Galerkin  were  done  on  a  Macllfx.  Each  run  takes 

approximately  30  min  and  1  Mbyte.  The  computations  for  the  GLS/ST  where  done  on  a  Cray 

■> 


4.  Conclusions 

We  have  compared  the  characteristic-Galerkin  and  Galerkin/ least-squares  space-time 
formulations  (which  are  both  suitable  for  time-dependent  domains!  based  on  error  estimates 
and  numerical  performance  for  a  test  problem  governed  by  the  advecuon-diffusion  equation. 
We  have  shown  that  the  formulation  of  the  problem  with  either  method  is  well-posed  and  the 
order  of  accuracy  is  the  same  as  when  the  spatial  domain  is  fixed,  that  is.  0{h  +  6r)  for  the 
characteristic-Galerkin  formulation  and  0(/i'  -!-5r’  ')  for  the  GLS/ST  formulation.  It  was 
brought  to  our  attention,  at  the  galley  proofs  time,  that  related  convergence  studies  appear  in 
[20,21],  The  test  problem  involves  the  transport  of  an  exponential  hill  in  a  rotational  flow 
field.  The  results  obtained  with  the  two  methods  for  various  combinations  of  the  Peclet 
number  and  mesh  motion  show  that  the  numerical  performance  of  these  two  methods  are  very 
good. 
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